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Statistical properties of the inhomogeneous version of the Olami-Feder-Christensen (OFC) model 
of earthquakes is investigated by numerical simulations. The spatial inhomogeneity is assumed to 
be dynamical. Critical features found in the original homogeneous OFC model, e.g., the Gutenberg- 
Richter law and the Omori law are often weakened or suppressed in the presence of inhomogeneity, 
whereas the characteristic features found in the original homogeneous OFC model, e.g., the near- 
periodic recurrence of large events and the asperity-like phenomena persist. 



Statistical properties of earthquake have attracted 
much interest in seismology as well as in statistical 
physics [HQ- It has been known for years that power- 
laws are often observed there, e.g., the Gutenberg- 
Richtcr(GR) law associated with the earthquake size dis- 
tribution or the Omori law associated with the time 
evolution of aftershock frequency, which leads to the 
view that earthquakes are essentially "critical" in nature 
P, H[ ■ A contrasting view of earthquakes has also been 
widespread, regarding earthquakes as "characteristic" in 
nature with characteristic energy and time scales [l|. 

In studying statistical properties of earthquakes, sim- 
plified models have been used, e.g., the so-called spring- 
block or the Burridge-Knopoff(BK) model The 
OFC model, which was first introduced by Olami, Feder 
and Christensen (OFC) as a further simplification of the 
BK model, is one of such statistical models of earthquakes 
11]. It is a two-dimensional coupled map lattice where 



the rupture is assumed to propagate from lattice site to 
its nearest-neighbor sites in a non-conservative manner, 
often causing multi-site seismic events or "avalanches" . 

In the past, many numerical studies have been per- 
formed on the OFC model, revealing that the model 
exhibits certain critical properties such as the GR law 
or the Omori law [171. More recent studies also 



unraveled the characteristic features of the OFC model 
3" 2(J. By investigating the time series of events, Ramos 



et al found the nearly periodic recurrence of large events 
[l8j . Kotani et al studied the spatiotemporal correlations 
of the model and identified in the OFC model a phe- 
nomenon resembling the "asperity" [I9j. Physical mech- 
anism underlying such asperity-like phenomena, i.e., a 
self-organization process of the stress concentration, was 
revealed in Ref.[20J. Thus, the OFC model, though an 
extremely simplified model, exhibits quite a rich phe- 
nomenology containing both critical and characteristic 
features observed in real seismicity. 

In the OFC model, "stress" variable /$ (/j > 0) is as- 
signed to each site on a square lattice with L x L sites. 
Initially, a random value in the interval [0,1] is assigned 
to each fi, while fi is increased with a constant rate uni- 
formly over the lattice until, at a certain site i, the fi 



value reaches a threshold, f c = 1. Then, the site i "top- 
ples", and a fraction of stress afi (0 < a < 0.25) is 
transmitted to its four nearest neighbors, while fi itself 
is reset to zero. If the stress of some of the neighboring 
sites j exceeds the threshold, i.e., fj > f c = 1, the site 
j also topples, distributing a fraction of stress afj to its 
four nearest neighbors. Such a sequence of topplings con- 
tinues until the stress of all sites on the lattice becomes 
smaller than the threshold f c . A sequence of toppling 
events, which is assumed to occur instantaneously, cor- 
responds to one seismic event or an avalanche. After an 
avalanche, the system goes into an interseismic period 
where uniform loading of / is resumed, until some of the 
sites reaches the threshold and the next avalanche starts. 
The transmission parameter a measures the extent of 
non-conservation of the model. The system is conserva- 
tive for a — 0.25, and is non-conservative for a < 0.25. 
A unit of of time is taken to be the time required to load 
/ from zero to unity. 

It should be noticed that the original OFC model is a 
spatially homogeneous model, where homogeneity of an 
earthquake fault is implicitly assumed. Needles to say, 
real earthquake fault is more or less spatially inhomoge- 
neous, which might play an important role in real seis- 
micity. Then, a natural next step might be to extend the 
original homogeneous OFC model to the inhomogeneous 
one where the evolution rule and/or the model parame- 
ters are taken to be random from site to site. 

Possible temporal variation of such a spatial inhomo- 
geneity might also be an important factor. There may 
be two distinct processes that could change the mate- 
rial parameters characterizing an earthquake fault from 
one seismic event to the next. The fast dynamical pro- 
cess during an earthquake rupture could change the fault 
state via, e.g., wear, frictional heating, melting, etc. In 
addition, in a long interseismic period until the next 
earthquake, the fault is subject to many slower processes, 
e.g., water migration, plastic deformation, chemical reac- 
tions, etc, which would necessarily cause the change in 
the material parameters at each position on the fault [l[ . 

In introducing the spatial inhomogeneity into the OFC 
model, there might be two extreme ways. In one, one may 
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assume that the randomness is quenched in time, namely, 
spatial inhomogeneity is fixed over many earthquake re- 
currences. In the other extreme, spatial inhomogeneity 
is assumed to vary with time in an uncorrelated way over 
earthquake recurrences. In the latter extreme, we assign 
spatial inhomogeneity independently to each successive 
event, resetting the one of the previous event. In order 
to get an overview of the role of spatial inhomogeneity in 
seismicity, full understanding of these two limits of the 
OFC model would be important and useful. 

There were several previous studies on the inhomo- 
geneous OFC model for both types of inhomogeneities 
3 21 "26 1- For the first type of inhomogeneity, i.e., the 
quenched or static randomness, Janosi and Kertesz intro- 
duced spatial inhomogeneity into the stress threshold and 
found that the inhomogeneity destroyed the SOC feature 
of the model [jl| . Torvund and Froyland studied the ef- 
fect of spatial inhomogeneity in the stress threshold, and 
observed that the inhomogeneity induced a periodic rep- 
etition of system-size avalanches [22|. Ceva introduced 
defects associated with the transmission parameter a, 
and observed that the SOC feature was robust against 
small number of defects [23| . Mousseau 24 1 and Bach et 



al 25| introduced inhomogeneity into the transmission 
parameter at each site. These authors observed that the 
bulk sites fully synchronized in the form of a system- wide 
avalanche over a wide parameter range of the model. 

For the second type of inhomogeneity, i.e., the dynam- 
ical randomness, Ramos considered the randomness as- 
sociated with the stress threshold, and observed that the 
nearly periodic recurrence of large events persisted [l8| . 
More recently, Jagla studied the same stress-threshold 
inhomogeneity, to find that the GR law was washed out 
by the small amount of randomness [26j . 

In the present paper, we study this second type of in- 
homogeneity with particular interest in the fate of the 
asperity-like phenomena identified in the homogeneous 
OFC model 0, [20[ • As to the form of the inhomogene- 
ity, there could be many different implementations. In 
the present paper, we consider the following three differ- 
ent forms of the spatial inhomogeneity and study their 
statistical properties by means of numerical simulations. 

Model[A]: Random anisotropy in the stress transmission 
When the site i topples, a stress ai/i is transmit- 
ted to its two nearest-neighbor sites and a stress o^/i is 
transmitted to the remaining two nearest-neighbor sites. 
Which neighbor is chosen to be a\ or a2 is decided at 
random at each site. 

Model [B]: Random magnitude in the stress transmission 
The transmitted stress is now isotropic among the four 
nearest neighbors of a toppled site, but its magnitude 
varies randomly from site to site. We assume that there 
are two possible values of the magnitude of the stress 
transmission parameters, a+ and a_. Which a being 



chosen is determined randomly at each site with the prob- 
abilities p and 1 — p, respectively. 

Model [C]: Random residual stress value 

When the site i topples, the stress fi is reset not to 
zero, but to a finite residual value of f r (0 < / r < 6 < 1), 
where f r is chosen uniformly between [0,6]. One may also 
consider the model where the threshold stress f c is not 
unity but distributes around unity (model[C']). Our pre- 
liminary study indicates that property of the model [C] 
is very much similar to that of model[C]. Hence, in the 
present paper, we concentrate on the model[C]. 

In our simulation, the lattice studied is L — 256 with 
open boundary conditions, and the pseudo-sequential up- 
dating proposed by Pinho et al being utilized [13]. Ini- 
tial 1 x 10 s avalanches are discarded to reach the steady 
states, and the subsequent 1 x 10 9 avalanches are gener- 
ated to study statistical properties. 

We begin with the properties of the model [A], i.e, 
the inhomogeneous OFC model where the transmission 
parameter exhibits site-random directional anisotropy, 
taking a value ai for two randomly chosen directions 
and a value a 2 for the remaining two directions. This 
randomness can also be characterized by its mean a = 
(ai + OL2)j2 and its standard deviation a = (a\ — a.2)j2. 

We show in Fig.l the size distribution of avalanches of 
the model [A] on a log-log plot for several values of the 
standard deviation a of the transmission parameter with 
a fixed mean a = 0.2. The avalanche size s is defined by 
the total number of "topples" in a given avalanche. As a 
is varied from (homogeneous model) to 0.015, the SOC 
feature of the original homogeneous model tends to be 
weakened, though a near-critical feature still persists to 
certain extent in the parameter range studied. 

The original homogeneous OFC model exhibits an- 
other well-known power-law, the Omori law (the inverse 
Omori law) associated with the time evolution of the fre- 
quency of aftershocks (foreshocks). A striking contrast 
to the homogeneous model occurs here. Fig. 2 exhibits 
on a log-log plot the time dependence of the frequency 
of aftershocks and foreshocks associated with mainshocks 
of s > s c — 100 for the case of a = 0.20 and a = 0.01. 
All events of arbitrary size which occur at an arbitrary 
site on the lattice just after (before) the mainshock of 
s > s c = 100 (t < 10~ 4 ) are counted as aftershocks 
(foreshocks) here. As can be seen from the figure, the 
introduced inhomogeneity destroys the Omori or the in- 
verse Omori laws almost completely, though a distinct 
power-law behavior was observed for the ho mog eneous 
OFC modelin the same i-range of Fig. 2 [13, US]- Here 
an increase or decrease of foreshocks/aftershocks itself is 
washed out by the introduced inhomogeneity. 

Next, we investigate the local recurrence-time distribu- 
tion P(T), which turned out to be an efficient probe of 
the characteristic feature of the model. The local recur- 
rence time T is defined as the time passed until the next 
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FIG. 1: (Color online) Log-log plot of the size distribution 
of seismic events of the model [A] for various values of the 
standard deviation a of the transmission parameter with a 
fixed mean value a — 0.20. The a = case corresponds to 
the homogeneous model. 
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FIG. 2: (Color online) The time dependence of the frequency 
of aftershocks and foreshocks of the model [A] on a log-log 
plot. The standard deviation of the transmission parameter 
is a = 0.01 with its mean a = 0.20. Mainshocks are the events 
of their size greater than s > s c — 100. No range threshold 
r c is imposed here in defining foreshocks or aftershocks. The 
time t is measured with the occurrence of a mainshock as an 
origin. 



avalanche occurs with its epicenter lying in a vicinity of 
the preceding avalanche, say, within distance r < r c (in 
units of lattice spacing) of the epicenter of the preceding 
event. In the original homogeneous OFC model, P(T) 
exhibited a sharp S- function-like peak at T = 1 — 4a [19j . 
The computed P(T) of the model[A] is shown in Fig. 3 
for various a and fixed a = 0.20. The size and the range 
thresholds are taken to be s > s c — 100 and r < r c = 10. 
As can be see from the figure, P{T) exhibits a clear peak 



structure even in the presence of inhomogeneity, indicat- 
ing that many events of the random OFC model[A] re- 
peat near periodically (the multi-peak structure of P(T) 
is originated from the size-threshold s c effect |2(|). As 
the c-value increases, the peak position of P(T) moves 
to a longer T-value, and the width of peak gradually in- 
creases: See the inset. 

In Fig. 4, we plot the position of the main peak, T*, 
of P(T) versus a for the cases of a — 0.19 and 0.20. As 
can be seen from the figure, T* increases with increasing 
a in proportion to a. In fact, the observed T*-value 
can be well fitted to the form T* = (1 - 4a)(l + 5a), 
as shown in Fig. 4. For the homogeneous case a = 0, 
this expression reduces to the one previousl y rep orted 
for the homogeneous model, T* = 1 - 4a [19|,l20j. The 
extra factor 1 + 5<r, which comes from the randomness, 
is closely related to the stress state of the model, as we 
shall see below. 
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FIG. 3: (Color online) Log-log plot of the local recurrence- 
time distribution of the model [A] for large avalanches of their 
size s > s c — 100. The standard deviation of the transmission 
parameter a is varied from (homogeneous model) to 0.15 
with a fixed mean- value of a = 0.2. The range parameter is 
r c — 10. The inset is a magnified view of the main peak. Note 
the difference in time scale from that of Fig. 2. 

In the original homogeneous OFC model, the sharp 
structure of P(T) arises due to the near-periodic rupture 
of the "asperity-like" cluster realized in the model 
[2pj , where nearly the same cluster of sites rupture many 
times near periodically with a period close to T* = 1 —4a 
with almost the same site as its epicenter. We have found 
that similar asperity-like events occur also in the present 
inhomogeneous model and are responsible for the sharp 
peak of P{T). An epicenter site of these asperity- like 
events lies close to the epicenter site of the preceding 
asperity event, but often moves from the previous one by 
several lattice spacings. This might be contrasted to the 
behavior in the homogeneous model where an epicenter 
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FIG. 4: (Color online) The peak position T* of the local 
recurrence-time distribution of the model [A] is plotted versus 
the standard deviation a of the transmission parameter for the 
mean-values of a — 0.19 and 0.20. The straight line represents 
the relation T* = (1 - 4a) (1 + 5a). 



site rarely moves during the asperity sequence 19|, . In 
the homogeneous OFC model, an epicenter site is located 
at the tip of the rupture zone where only one out of four 
nearest-neighbor sites is contained in the rupture zone 
(rij = 1 site), while it is never located at the interior 
site inside the rupture zone (n.j = 4 site) [2(J. In the 
present inhomogeneous model, an epicenter site tends to 
be located at the corner of the rupture zone where two out 
of four nearest-neighbor sites are contained in the rupture 
zone (jii = 2 site) , while a finite fraction of epicenter sites 
are located at the interior site inside the rupture zone. 

In Fig. 5, we show the stress distribution D(f) at the 
time of toppling for the sites in the asperity cluster. D(f) 
exhibits a peak at a stress value f = f p exceeding the 
threshold f c = 1, whose position depends on the a- value. 
In the inset of Fig. 5, we plot the observed f p - values as 
a function of a. The data are well fitted to the form 
f p = 1 + 5<7, as shown in the inset. The observation 
that the asperity site topples at a stress value around 
fp = 1 + 5c in the inhomogeneous model, rather than at 
/ = 1 as in the homogeneous model, might explain the 
observed deviation of the recurrence time of the asperity 
events from that of the homogeneous model. Naturally, 
the random modulation in the transmission parameter 
a of order a would suppress the formation of the stress 
state highly concentrated at the threshold / = 1 by an 
amount of a. Meanwhile, the reason why this factor is 
given by 1 + 5cr quantitatively is not clear to us. 

Next, we study the model [B], i.e., the inhomogeneous 
OFC model where the transmission parameter, though 
isotropic in its direction, exhibits a random distribution 
in its magnitude from site to site, taking a value a+ with 
the probability p and a value a_ with the probability 
1 — p. This randomness can also be characterized by 
its mean value, a = pa + + (1 — p)a-, and its standard 
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FIG. 5: (Color online) The stress distribution D(f) of the 
model [A] at the time of toppling of each site contained in 
the rupture zone of the asperity events for several values of 
the standard deviation a of the transmission parameter. The 
mean value is fixed to a = 0.2. The inset represents the 
peak position of D(f) plotted versus a, while the straight 
line represents the relation f p = 1 + 5a. 



deviation, a = \Jp{\ — p)(a+ — a-). We have studied 
various cases of aq_, a_ and p, and have found that the 
properties of the model [B] is determined by its mean a 
and standard deviation a. 

Indeed, the properties of the model [B] turns out to be 
very much similar to those of the model [A] . The SOC fea- 
ture of the size distribution tends to be weakened, while 
the foreshocks/aftershocks activity is suppressed almost 
completely. In particular, the local recurrence-time dis- 
tribution of the model [B] exhibits a clear peak structure 
borne by the asperity-like events at T = T* which are 
very well fittable by the relation T* = (1— 4a)(l+5ir), the 
same formula as we obtained for the model [A]. Hence, 
the recurrence time of the asperity events of the homoge- 
neous OFC model and the inhomogeneous OFC models 
[A] and [B] with the random transmission parameter is 
given by a common simple formula, T* = (1— 4a)(l+5a). 
As in the model [A], the stress distribution at the time 
of toppling for the sites in the asperity cluster also ex- 
hibits a peak at a stress value / = f p > fc = 1. Again, 
the observed / p -values are well fittable by the formula 
f p = 1 + 5er, explaining the origin of the observed recur- 
rence time of the asperity events, T* = (1 — 4a)(l + 5a). 

Finally, we study the model [C], i.e, the inhomoge- 
neous OFC model where the residual stress value exhibits 
a site-random distribution, taking uniformly a value be- 
tween [0,(5]. The transmission parameter a is assumed 
to be homogeneous here. Concerning the avalanche size 
distribution and the time evolution of the frequency of 
foreshocks/aftershocks, a tendency very much similar to 
the ones of the models [A] and [B] is observed, i.e.. the 
SOC feature of the size distribution tends to be weak- 
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ened, while the foreshocks/aftershocks activity is sup- 
pressed almost completely. 

Fig. 6 exhibits the local recurrence-time distribution of 
the model [C] for several values of S with fixed a = 0.20. 
P(T) exhibits a clear peak structure as in the cases of 
models [A] and [B] , whereas the main peak now lies pre- 
cisely at T* = 1— 4a irrespective of the value of S, though 
the peak is broadened with increasing 5. Thus, the ran- 
domness in the residual stress does not affect the recur- 
rence time of the asperity events, in contrast to the ran- 
domness in the transmission parameter. The peak events 
are borne by the asperity-like events. An epicenter site 
here tends to be located either at the corner (n* = 2) or 
at the boundary (rij = 3) of the rupture zone. 
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FIG. 6: (Color online) Log-log plot of the local recurrence- 
time distribution of the model [C] for large avalanches of their 
size s > s c = 100. The residual-stress threshold S is varied 
from (homogeneous model) to 0.1, while the transmission 
parameter is taken to be uniform a = 0.2. The range param- 
eter is r c — 10. The inset is a magnified view of the main 
peak. 

In Fig. 7(a), we show the stress distribution D{f) of 
the model [C] at the time of toppling for the sites in the 
asperity cluster. D(f) exhibits a peak at a stress value 
/ = f p > f c = 1, whose position depends on the 6- value. 
Since the residual stress is nonzero and site-dependent, 
the distribution of the stress drop A/ at each site shown 
in Fig. 7(b), rather than the one of the stress itself at the 
time of toppling shown in Fig. 7(a), might exhibit a peak 
at A/ = 1. Indeed, this seems to be the case as can 
be seen from the figure. This observation for the stress 
drop is fully consistent with our observation in Fig. 6 that 
the recurrence time of the model[C], T* = 1 — 4a, is 
independent of S. 

In summary, we studied the statistical properties of the 
inhomogeneous version of the Olami-Feder-Christensen 
(OFC) model by numerical simulations. The spatial in- 
homogeneity was assumed to be dynamical, varying from 
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FIG. 7: (Color online) The stress distribution D(f) of the 
model [C] at the time of toppling of each site contained in 
the rupture zone of the asperity events. The transmission 
parameter is a — 0.2, while the residual-stress threshold is 
either 5 — 0.05 or 0.1. In (a), the stress value at the time 
of the toppling is taken as the abscissa, whereas, in (b), the 
stress drop, i.e., the stress value at the time of the toppling 
minus the residual-stress value at that site, is taken as the 
abscissa. 



event to event in an uncorrclatcd way. In common with 
all three types dynamical inhomogeneities studied, the 
critical features of the original homogeneous OFC model 
are often weakened or suppressed: With increasing the 
inhomogeneity, the deviation from the GR law becomes 
more pronounced, accompanied by the suppression of 
large events. Our result corroborates the recent obser- 
vation of Jagla for the model [C], a model similar to our 
model [C] [26( . (We note in passing that Jagla further re- 
ported that the inclusion of the slow structural-relaxation 
process in the inhomogeneous OFC model [C] apparently 
revived the GR law [26|.). The foreshock/aftershock ac- 
tivity described by the Omori (or inverse Omori) law is 
entirely gone. By contrast, the characteristic features of 
the original homogeneous OFC model persist: In all types 
of dynamical inhomogeneities studied, near-periodic re- 
currence of large events persist borne by the asperity- 
like events. We emphasize that characteristic features 
are observed for all inhomogeneities in common, suggest- 
ing that the persistence of the characteristic features is a 
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generic property of the dynamical inhomogeneity. 

Note also that the properties of the dynamically inho- 
mogeneous models are quite different from those of the 
static or quenched inhomogeneous models. In the latter 
case, the introduced inhomogeneity often gives rise to a 
full synchronization and a periodic repetition of system- 
size events. Such a system-wide synchronization is never 
realized in the present dynamically homogeneous mod- 
els. Presumably, temporal variation of the spatial inho- 
mogeneity may eventually average out the inhomogeneity 
over many earthquake recurrences, giving rise to the be- 
havior similar to that of the homogeneous model. 
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